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A THREE-DIMENSIONAL VISCOUS/POTENTIAL FLOW 
INTERACTION ANALYSIS METHOD FOR MULTI-ELEMENT WINGS 

By 

F.A. Dvorak, F.A. Woodward and B. Maskew 
Analytical Methods, Inc. 


SUMMARY 


An analysis method and computer program have been developed 
for the calculation of the viscosity-dependent aerodynamic charac- 
teristics of multi-element, finite wings in incompressible flow. 
This work is an extension to three dimensions of the method deve- 
loped previously under Contract NAS2-7048, Reference _1. The 
methods differ in that a fully three-dimensional potential flow 
program is now used to determine the inviscid pressure distribu- 
tion about the configuration. The potential flow program uses 
surface source and vortex singularities to represent the inviscid 
flow. The method is capable of analysing configurations having at 
most one slat, a main element, and two slotted flaps. Currently, 
the configurations are limited to full span slats or flaps. The 
configuration wake is allowed to relax as a force-free wake, al- 
though roll-up is not allowed at this time. 

Once the inviscid pressure distribution is calculated, a 
series of boundary layer computations are made along streamwise 
strips. Each strip is treated as if it were a separate infinite 
yawed wing; consequently, the boundary layer methods developed 
in Reference 1 are used directly in the new program. Source dis- 
tributions are determined along the streamwise strips from the 
boundary layer and potential flow calculations. They are used to 
modify the boundary conditions in the second and subsequent cal- 
culations of the potential flow. These sources represent the 
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effect of the boundary layer in the modification of the potential 
flow. The sequence of potential flow and boundary layer calcula- 
tions continue until convergence of the lift coefficient (usually 
2 to 4 iterations) . Lift, drag and pitching moment coefficients 
are then determined. The method is currently the only one of its 
kind capable of analyzing finite wings. 

The computer program is written in Fortran IV for the CDC 
6600 and 7600 family of computers. The program operates in the 
overlay mode, and requires an amount of storage dependent on the 
number of panels used to describe the configuration. Typically 
300 panels will require about 170,000 (Octal) words of storage. 
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INTRODUCTION 


Background 

An analysis method has recently been developed which is capa- 
ble of predicting the aerodynamic characteristics of infinite 
yawed multi-element wings (Ref. 1) . The effect of viscosity is 
included in the calculation through the use of distributed sources 
determined from the boundary layer analysis. These sources are 
included in the subsequent potential flow calculation giving a 
viscosity dependent flow field. Convergence of the method as 
measured by lift coefficient is generally very rapid. 

Application of the method for comparison with experiment has 
demonstrated its validity for a wide range of problems. One such 
case is the optimization of flap-gap and overlap for a multi-ele- 
ment configuration resulting in maximum lift for a given angle- 
of-attack. An extended version of the method (Ref. has been 
coupled with an optimization program for use in determining the 
optimum slat locations for maximum lift. 

Because most aerodynamic flows of interest are three dimen- 
sional in nature, the ability to consider the effect of finite 
span on such parameters as flap-gap is very desirable. A fully 
three-dimensional method would have many applications important 
to the des:.gner of modern high-lift systems. Its availability 
would greatly reduce the reliance on the experimental method. 
Development of such a method was consequently undertaken, and the 
resulting procedure is described in this report. 

, Problem Definition 

The calculation of the potential flow field about a general 
three-dimensional multi-element wing (see Figure 1) represents 
the first task of any analysis method. Because of the complex 
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nature of the configurations likely to be analysed, a geometry 
control program will be an essential feature of the overall method. 
As in the case of Program VIP (Ref. 1) , the potential flow method 
must be capable of predicting the pressure distribution at selec- 
ted off-body points above the flap elements. An additional requ- 
irement in the three-dimensional case is a wake representation 
capable of providing for wake roll-up and ultimately for tip-edge 
vortex roll-up. In three dimensions these effects become impor- 
tant at relatively low lift coefficients. 

With the potential flow field specified, it is necessary to 
predict the boundary layer development over the multi-element con- 
figuration. Calculations must include stagnation line initial 
conditions, laminar, transition and turbulent boundary layer deve- 
lopments and laminar or turbulent separation predictions for each 
element of the wing high- lift system. The calculations must in- 
clude accurate predictions of boundary layer development in the 
regions where wing or first flap upper surface and cove boundary 
layers merge with the downstream flap upper surface boundary layer. 
Both longitudinal curvature and normal pressure gradient terms 
must be included in the governing boundary layer equations as each 
effect has a significant influence on the boundary layer develop- 
ment, and subsequently on the section drag coefficient. These 
effects are particularly important in the wing trailing-edge/f lap 
leading-edge region. Once the boundary layer development is known, 
its effect on the external flow must be determined. 

A complete analysis program for the aerodynamic characteris- 
tics of general three-dimensional multi-element wings is developed 
by combining the, separate goemetry, potential flow and boundary 
layer calculation procedures. Iteration between the separate pro- 
grams results in the prediction of viscosity-dependent aerodyna- 
mic forces. The different calculation schemes that form the 
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elements of the integrated computer program are discussed in the 
following sections. 
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POTENTIAL FLOW METHOD 


Configuration Definition 

A typical panel subdivision of a high-lift wing configuration 
is illustrated in Figure 1. A reference coordinate system is esta- 
blished with origin at or near the leading edge of the configura- 
tion, having an x axis lying in the plane of symmetry parallel to 
the streamwise axis, and a vertical z axis. 

Each wing element may be specified in its own or in the refe- 
rence coordinate system. Individual coordinate systems are related 
to the reference coordinate system by pivot points. The pivot 
points are prescribed in the reference coordinate system. In order 
to loft the configuration, element rotation angles must also be 
prescribed. Given the pivot points and rotation angles, any ele- 
ment may be translated and rotated to the desired location relative 
to the reference coordinate system. 

Wing twist, taper and dihedral can be readily accommodated by 
the geometry control program. 

If the wing configuration is made up of a main element and 
one or more slotted flaps, then additional analysis is required 
to determine flap upper surface longitudinal radius or curvature 
for later use in the finite difference boundary layer calculations. 
The spline technique as described ir. Reference ^ is used in the 
present method. 

Each element of a multi-element wing configuration is repre- 
sented by a series of streamwise airfoil sections given at selec- 
ted intervals along the span. Planar surface panels are located 
between adjacent sections, with the comer points being defined 
by the section coordinates. The section coordinates may be given 
in percent chord or directly in terms of the reference coordinate 
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system. The panels in each wing section are generally numbered 
sequentially from the trailing edge on the lower surface around 
the leading edge to tlie trailing edge on the upper surface. The 
same number of points at approximately the same percent chord 
locations must be used to define the wing upper and lov;er surfaces 
Each chordwise strip of wing panels contains a constant vor- 
tex distribution starting at the trailing edge of the lower sur- 
face, around the leading-edge, and ending at the trailing edge of 
the upper surface. These vortex strips are used to provide circu- 
lation around lifting srurfaces. One additional panel is defined 
in the wake aft of each chordwise strip to provide control points 
for satisfying the Kutta condition. The additional panel lies in 
the plane of the trailing edge bisector. 


z 
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Each panel is defined by foui* corner points. Since these four 
points may not lie in the same plane, an equivalent planar panel 
is generated using the method of Reference _3, the details of v?hich 
are described in Appendix I. A panel coordinate system is defined 
which has its origin at the panel centroid. The S and n axes lie 
in the plane of the panel, v;hile the C axis is perpendicular to 
that plane. The C axis is oriented such that the C plane is 
parallel to the reference x axis. 

C 



Figure 2 . Panel Coordinate System 


Since the velocity components induced by the source and vor- 
tex distributions are given in terms of the panel coordinate sys- 
tem, a nine-element transformation matrix, ^ ^ j / is calculated for 

each panel to transform the coordinates of points and the compo- 
nents of vectors from the reference coordinate system to the panel 



coordinate system. In addition, the panel area, the coordinates 
of the centroid, and the length of the principal diagonal are cal- 
culated. 


Inviscid Flow Model 

The surface of the configuration is subdivided into a large 
number of panels, each of which contains a constant source and 
constant vortex distribution. A concentrated vortex originates 
from the downstream corners of each vortex panel, and follows the 
wing contour along the inner and outer edges of each strip of pa- 
nels until the trailing-edge is reached. Downstream of the trail- 
ing-edge, this system of vortices forms the trailing vortex wake. 

Analytical expressions for the perturbation velocity field 
induced by a constant source distribution on an arbitrary quadri- 
lateral panel are given by Hess and Smith (Ref. 3 ) . Equivalent 
expressions for the perturbation velocities induced by a constant 
vortex distribution on the same arbitrary quadrilateral are given 
herein. The perturbation velocities are used to calculate the 
coefficients of a system of linear equations relating the magnitude 
of the normal velocities at the panel control points to the un- 
known source and vortex strengths. The source and vortex strengths 
which satisfy the boundary condition of tangential flow at the 
control points for a given Mach number and angle of attack are 
determined by solving this system of equations by an iterative 
procedure. Initially, the vortex wake is constrained to trail in 
the free stream direction during this solution procedure. However, 
an optional wake relaxation procedure is also available which 
allows the trailing vortices to move in a vertical direction and 
approximate the streamlines of the flow. Transverse movement of 
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the vortices is not permitted in order to avoid adverse interfe- 
rence between the vortices and control points on downstream flap 
elements. With this option, an additional iteration loop is super- 
imposed on the basic solution procedure, in which the influence of 
the wake on the wing is recalculated during each cycle. Complete 
wake relaxation, including roll-up of the wing tip vortices, is 
not permitted in order to minimize the number of iterations re- 
quired for convergence. 

The pressure coefficients at panel control points are calcu- 
lated from the perturbation velocity components, and finally, the 
forces and moments acting on the complete configuration are ob- 
tained by numerical integration. 


The Incompressible Velocity Components 

The perturbation velocity components, u , v , and w , induced 
by a constant source distribution on an arbitrary quadrilateral 
panel are derived in Reference _3. 

The equations are given in terms of the panel coordinate 



Figure 3. Panel Geometry 
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The panel corner points are numbered in a clockwise direc- 
tion. The perturbation velocities at an arbitrary point, P(x,y,z), , 
are given as the sura of the contributions of the four sides of the 
quadrilateral as follows. 
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- 1 

F. = tan - 
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Clearly, the expressions for the velocity components given above 
may be extended to a panel having any mjmber of sides. 

For points located at a distance from the centroid greater 
than four times the length of the major diagonal, the quadrilate- 
ral is approximated by a point source at the centroid. In this 
case, the .expressions for the velocity components are considerably 
simplified, becoming 

u = 2A ’(4) 

® d3 


V 


S 



(5) 


w 


s 


zA 

d3 


( 6 ) 


where 


2 2 2 
= X + y + z 


and 


A = panel area 
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Additional multiple expansions given in Ref. ^ for the velocity 
components of points located at intermediate distances from the 
centroid are not used in this program. 

The perturbation velocity components, u^, v^ and w^, induced 

by a constant vortex distribution on the same arbitrary quadrila- 
teral panel may be derived in a similar manner. In this case, the 
direction of the elementary bound vorticity must be specified in 
advance, and is chosen to be parallel to the panel leading edge 
(between corners 1 and 2) • 


TRAILING 

VORTEX 

SHEET 



Figure 4. Vorticity Distribution 

The contribution of the bound vorticity is obtained by integration 
(see Ref. 4^) . The resulting perturbation velocities at an arbi- 
trary point, P(x,y,z) are obtained: 
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( 7 ) 


u = w 
V s 


V„ = -B 
V i 2 s 


w = B V - u 
V 12 s S 


( 8 ) 


(9) 


The bound vorticity on the panel generates linearly varying 
vortices along the side edges 1-4 and 2-3. These edge vorti- 
ces have zero strength at the leading edge , and a strength at 

the trailing edge, where ^ The 

edge vortices continue downstream from the trailing edge, main- 
taining a constant strength, ^ . If the side edges are of un- 
equal length, a constant vortex sheet is contained between the two 
trailing vortices having a net strength equal to their difference 
in vorticity. 

Each of these edge and trailing vortices makes a contribution 
to the perturbation velocities. Formulas for these velocity in- 
crements are derived in Ref. 4., and tabulated in Appendix II. 


Compressibility Corrections 

The velocity components in compressible flow are found by 
applying Gothert's Rule (Ref. 5). Two options are available in 
the program for applying the compressibility corrections, and are 
designated Rule 1 and Rule 2. 

Rule 1 applies the method originally proposed by Gothert. 

The incompressible velocity components are calculated on an analo- 
gous body obtained by the following transformation; 
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X 


= By (10) 

z = Bz 

a 

where 

B = 

The boundary conditions of tangential flow are applied on the an- 
alogous body, and the resulting incompressible perturbation velo- 
cities are transformed back to the real body by 

u = u,/B^ 

ci 

V = v^/B (11) 

w = w^/B 

a 

The total velocity vector at a given point is then 
= Uflo cos a cos 3 + u 

V = sin 3 + v (12) 

W = U sin a cos 3 + w 

00 

It is now known that this compressibility rule yields good results 
only for slender bodies at small angles of attack. The validity 
of this rule decreases with increasing values of the surface slope. 
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This effect is particularly noticeable for tv70-dimensional airfoil 
sections. In the vicinity of the nose, Gothert ' s Rule (which is 
equivalent to the Prandtl-Glauert Rule in this example) gives ex- 
cessively high suction peaks on the upper surface. The reason for 
this failure of the theory is the manner in which the boundary con- 
ditions are satisfied at the surface of the analogous body which 
is thinner by the factor B than the real body, the curvature of 
the flow near the nose is correspondingly increased, resulting in 
higher suction peaks. In order to eliminate this effect, it is 
necessary to satisfy the boundary conditions on the surface of the 
real body. 

Rule 2 was first proposed by Kraus in Reference 6. Beginning 
with the analog body as before, the expressions for the perturba- 
tion velocity components are corrected for compressibility, using 
Equation (11), prior to solving the boundary condition equations. 
The boundary conditions of tangential flow are then applied on the 
surface of the real body, resulting in improved results for the 
velocities and pressure coefficients. 

The Boundary Condition Equations 

The boundary condition of tangential flow at panel control 
points establishes a system of linear equations for determining 
the strengths of the source and vortex distributions. The geo- 
metrical relationship between each panel and control point is re- 
quired to evaluate the coefficients of this system of equations 
for a given free-stream Mach number. 
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Normal Velocity at Panel Control Points 

Each surface panel is assigned a control point located at the 
panel centroid. Each vortex strip is assigned a control point just 
behind the trailing edge of the wing in the plane of the trailing 
edge bisector. (This point is normally located 0.1 percent of the 
local chord behind the trailing edge.) 

The resultant velocity normal to panel i at its control point 
the sum of the normal component of the free-stream velocity 
vector and the normal velocities induced by the panel source and 
vortex distributions. Setting the magnitude of the free-stream 
velocity vector equal to unity, its component normal to panel i is 

R. = cos a*cos 0*n + sin g*n + sin a‘cos B*n (13) 

1 X. y. z. 


where n , n , and n are the direction cosines of the normal of 

X • y 4 z « 

panel i (see Appendix I ) , a is the angle of attack and 3 is the 
angle of yaw of the free-stream velocity vector in the reference 
axis system. 

The normal component of velocity induced at the control point 
of panel i by the source and vortex distributions is given by 




+ 





(14) 


where v , v , and v are the three components of velocity 
X . . v . . z . . 

13 ^13 13 

parallel to the reference axis at control point i induced by a 
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unit strength source or vortex distribution on panel j , is the 

aerodynamic influence coefficient, and Oj is the strength of the 
singularity. 

The three components of velocity parallel to the reference 
axes are obtained by multiplying the velocity components calculated 
in the panel coordinate system by the transformation matrix given 
in Appendix I. For example, 
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^ij 

= u . . 
13 
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13 

+ 
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13 
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w. . 
13 
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13 
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13 

n 

z . . 
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(15) 


Combining Equations (13) and (14) , 


n . 

1 


R. + V 


n 


R. + 
1 . 


Z a . .a . 
13 3 

j=l 


(16) 


Solution of the Boundary Condition Equations 

The boundary condition of tangential flow at panel control 
points is satisfied if the normal velocities are set equal to zero 
on all panels. 
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Thus 


0 


i = 1,N 


HijOj = -R^ . i = 1,N (17) 

In matrix notation, 

[A^j] {a^} = -{ R^} (18) 

where A^j is the matrix of aerodynamic influence coefficients, and 

the right side of the equation is given by Equation (13). 

This sytem of equations can be solved by direct inversion to 
determine the unknown source and vortex strengths. However, for 
the large order matrices usually encountered in aerodynamic prob- 
lems, an iterative solution procedure described in Reference 1 _ is 
used. 

The aerodynamic matrix is subdivided into smaller blocks, 
with no block exceeding order 60. The matrix elements in each 
block are carefully chosen to represent some well-defined feature 
of the configuration. For example, a wing block represents the 
influence of one chordwise column of wing panels, and includes 
the influence of the vortex strip as well as the source panels. 

The initial iteration calculates the source and vortex 
strengths corresponding to each block in isolation. For this step^ 
only the diagonal blocks are present in the aerodynamic matrix, 
and the solution is obtained by a direct inversion of the diagonal 
blocks. Once the initial approximation to the source and vortex 


or 


ni 


N 

Z 

j=l 
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strengths is determined, the interference of each block on all the 
others is obtained by matrix multiplication. The incremental nor- 
mal velocities obtained are subtracted from those specified on the 
right side of Equation (18) . This process is repeated until the 
residual interference velocities are small enough to ensure that 
convergence has been achieved. 

In the method of Reference four optional iteration proce- 
dures are available to provide rapid convergence. These are: 

1. Blocked Jacobi (JB) 

2. Blocked Gauss-Seidel (GSB) 

3. Blocked Successive Over- relaxation (SORB) 

4. Blocked Controlled Successive Over- relaxation (CSORB) 

In general, CSORB gives the most rapid convergence. 


Calculation of the Pressures, Forces and Moments 

Once the source and vortex strengths have been determined, 
the three components of velocity at control point i may be obtained 
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cos a cos 8 + 



(19) 


V. 

1 
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i 


sin B + 


sin a cos B + 



( 20 ) 


( 21 ) 


where the includes both source and vortex strengths, and 
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are defined following Equation (14) . The pres- 



The forces and moments acting on the configuration can now be 
obtained by numerical integration. The normal force, side force, 
axial force, and pitching moments (about the origin of coordinates) 
of panel i are given by: 
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(28) 


_ V -7 - Z.X; 

^ i 


= Y,., - Vi 


(29) 


and n are the direc- 

here A is the area of the panel, z. 

' t andx y. ana z, axe the coorainatea 

tion cosines of the normal, and x., 

of the panel control point. ^„^„i^ients are obtained by snm- 

total force and of the plane of 

Ohe panel forces and moments on both 

symmetry : 
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( 35 ) 



A c 
w 





Finally, the lift, side force and drag coefficients are: 

cos a - (C^ cos B “ 

Cg = cos 6 + sin B (37) 

Cjj = (Cjj cos 3 - sin B) cos a + sin a (38) 

The program computes the forces and moments acting on the 
wing, and sums them to obtain the total forces and moments on the 
complete configuration. Wing section forces and moments may also 
be calculated at the user's option. 
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Relaxed Wake Model 


In the analysis of configurations with high circulation, 
the orientation of the trailing vortex sheet can have an important 
effect on the calculated results, particularly in the case of mul- 
tiple components, e.g., a wing with slotted flap or wing with 
tailplane. The trailing vortex sheet should, in fact, carry no 
load, but, satisfying this condition renders the problem non-linear 
insofar as the pressure distribution and the wake shape are inter- 
dependent. Analysis of such problems, therefore, requires an 
iterative procedure in which the wake shape is first assumed in 
order to solve for the pressure distribution; knowledge of the 
pressure distribution allows the streamlines to be calculated, 
which give the force-free orientation of the trailing vortex sheet. 
Further iterations may be necessary to make the pressure distri- 
bution and wake shape compatible. 

In the present method, the wake model and iterative tech- 
nique are similar to those used in a quadrilateral vortex- lattice 
method (Ref. £ ) . The trailing vortex sheet is represented by a 
system of discrete vortices (Figure 5 ) attached to the panel edges 
at the trailing-edge of each component (i.e., wing, flap, etc.). 
Each vortex is segmented over the region from its starting point 
to a station a short distance downstream of the last component. 

The segment lengths in this region can be varied from vortex to 
vortex so that small segments can be used where the vortices are 
expected to bend most, and larger segments can be used where little 
displacement is expected. This helps to minimize computing time. 
Downstream of the segmented region, each vortex is straight, semi- 
infinite and streamwise. Initially, all the vortices are assumed 
straight in the streamwise direction. In the iterative procedure, 
each vortex segment is aligned with the local mean velocity — 
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CIRCULATION VALUES ON THE LEFT 
AND RIGHT CHORDWISE STRIPS 



MEAN VELOCITY GIVES 
DIRECTION FOR NEXT SEGMENT 


Figure 5. Trailing Vortex Segments 
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calculated near the segment mid point. Reference Figure 5 - to 
render the segment approximately force free. But, the spanwise 
movement of the vortices is not allowed in the present program; 
this restriction is necessary at this time to avoid spurious re- 
sults when the discrete trailing vortices from one component pass 
close to the surface control points on another surface. Canceling 
the spanwise movement of the vortices maintains the correct rela- 
tionships between the discrete vortices and the control points. 
This is a near-field problem associated mainly with discrete vor- 
tex models, and is particularly apparent in clcse-approach situa- 
tions such as when the wing wake passes over a slotted flap. How- 
ever, techniques such as those of Refs. 9 and 10 can remove this 
restriction to allow a complete wake roll-up calculation. 

The strength of each trailing vortex is the difference in 
circulation between the two chordwise panel strips adjacent to the 
start of the vortex (Figure 5 ) . The strengths are changed at 
each iteration (Figure 6 ) . To minimize the computation in each 
iteration, the influence coefficient matrices for the wake vor- 
tices are stored separately. Clearly, the surface panel influence 
coefficients remain constant as the wake shape changes; only the 
wake contributions need be recalculated in each iteration. Stor- 
ing the wake influence coefficients separately allows the new set 
to be compared with the previous set; then only the differences 
are added to the main influence coefficient matrices. 

The number of iterations for wake shape is an input parameter 
at this stage. The first iteration gives most of the non-linear 
effect, and is probably as far as the calculation need be taken 
with the present restriction of canceling the spanwise movement. 
This model should be adequate up to values of the order of 2. 

Higher lift values, and also more detailed evaluation of pres- 
sures such as in the tip-edge vortex region, will require the 
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Figure 6. Wake Iteration Procedure 








complete roll-up calculation and more iterations. In this case, 
a convergence criterion would be required to control the number 
of iterations. This criterion could be based on the changes in 
wake influence coefficients or on some overall parameter such 
as Cj^. 
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BOUNDARY LAYER CALCULATION METHODS 


The finite span multi-element wing is divided into a number 
of streamwise strips. Each strip is treated as if it were a 
separate infinite span wing. On each strip the boundary layer 
development is calculated along stream lines from the stagnation 
line to the trailing edge of each element using the boundary 
layer methods currently in Program VIP (Ref. . No major modi- 
fication of the boundary layer programs has been made. At the 
current state of development of three-dimensional analyses for 
high-lift systems, it is believed that the infinite span boundary 
layer methods adequately predict the boundary layer development. 

At high angles-of-attack, or with very low aspect ratios the 
infinite-span methods are not expected to accurately predict the 
boundary layer development, particularly in the wing tip region. 
However, other phenomena, such as tip vortex formation, wake roll- 
up and main wing vortex wake-flap interaction greatly influence 
the inviscid flow, and must be taken into account before fully 
three-dimensional boundary layer methods are considered. 

The boundary layer methods are described in detail in Refe- 
rence and are summarized briefly in the following paragraphs . 

Stagnation Line Flow 

Theoretical predictions of the stagnation line flow of an 
infinite yawed wing by Cumpsty and Head (Ref. and Bradshaw 

(Ref. as well as others indicate that the boundary layer 

approaches an asymptotic state where frictional forces are balan- 
ced by divergence of the flow from the spanwise to the streamwise 
direction. Cumpsty and Head found that the stagnation line bound- 
ary layer integral parameters (H, 9, and C^) and the state (laminar 
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or turbulent) correlate with the parameter C* = V^/ (vdU) . 

ds 

Cumpsty and Head (Ref. 3^) later experimentally verified their 
theoretical correlations, and it is these correlations that are 
used to determine the boundary layer characteristics on each 
element of the swept finite wing.. If the wing is unswept, then 
conventional two-dimensional correlations are used to start the 
boundary layer calculations. 

Conventional Boundary Layer Methods 

Integral boundary layer m.ethods are used for all conventional 
boundary layers such as on the upper and lower surfaces of the 
main element, the lower surfaces of flap elements and the upper 
surfaces of flap elements up to the slot exits. If the Reynolds 
number is sufficiently low to allow laminar flow on a swept wing, 
the two-dimensional equations of Curie (Ref. 1_4 ) are solved along 
external streamlines to determine the laminar boundary layer 
development. It is assimied that laminar cross-flow effects have 
a negligible influence on the overall calculation, at least for 
moderate sweep angles. 

The streamwise boundary layer characteristics are used with 
the correlation of Smith (Ref. 1^) to determine the point of la- 
minar instability. With the point of instability, the momentum 
thickness Reynolds number distribution, Rg, and the pressure gra- 
dient parameter, k, known, the transition point is determined 
using Granville's correlation (Ref. ^). The turbulent boundary 
layer development over an infinite swept wing is calculated using 
the method of Cumpsty and Head (Ref. 1^). If the initial stagna- 
tion line flow is turbulent, Ciampsty and Head's method is used 
from the stagnation line to the element trailing edge. 
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In those cases where laminar separation occurs prior to 
transition, a correlation based on the data of Caster (Ref. 18 ) 
is used to determine if turbulent reattachunent occurs. Should 
reattachment be predicted, the calculation continues for turbulent 
flow; if catastrophic separation is predicted, the boundary layer 
calculation is terminated. 

Confluent Boundary Layer Method 

The boundary layer development from the slot exit to the 
trailing edge of each flap upper surface is determined using an 
infinite swept wing version of a finite difference method deve- 
loped by Dvorak (Ref. This method treats the upstream eleitmt 

wake and flap boundary layer as if it were all one thick boundary 
layer with initially an embedded potential core. The differencing 
scheme is based on the method of Crank and Nicholson (Ref. 20 ) , 
as modified by Dvorak and Head (Ref. ^). The turbulent closure 
is via an eddy viscosity model modified for the confluent bound- 
ary layer from the two-dimensional model for bo\indary layers and 
wall jets developed by Dvorak (Ref. 1^). The calculations include 
the effects of longitudinal surface curvature and the variation 
of static pressure through the boundary layer. The static pres- 
sure field, p(s,z), is determined directly from the potential 
flow solution. 

The initial velocity profile required to start the finite 
difference calculation at the slot exit is constructed from: (1) 

the integral boundary layer solution at the slot exit on the 
upper surface of the component in question; (2) the potential 
core as determined from the potential flow solution; and (3) the 
upper and lower surface boundary layer solutions at the trailing 
edge of the upstream element. With the initial velocity profile 
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the surface curvature and the static pressure field known, the 
boundary layer equations are solved in a forward inarching fashion 
to the trailing-edge of the component. 


VISCOUS/INVISCID INTERACTION 


The effect. of boundary layer displacement on the potential 
flow is simulated by distributing sources of known strength on 
the panels used to describe the wing geometry. The strengths of 
these sources are determined directly from the boundary layer 
solutions as q. = - 3 — (U. 6 .) where U. is the streamwise potential 

flow velocity at the edge of the boundary layer and 6 ^ is the 

streamwise displacement thickness. 

The addition of this source distribution modifies the normal 
velocity at the control point of panel i. Consequently, the bound 
ari' condition (Equation 16) is modified as follows 

N 

Vni = 

j = l 

Since q^ is known for each iteration, the right hand side only of 
Equation 18 is altered, giving (in matrix notation) 



(40) 


Because the original geometry is not modified by the use of 
distributed sources, the aerodynamic influence coefficient matrix 
need not be recalculated. Subsequent iterations between the po- 
tential flow and boundary layer calculations results in convergent 
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solutions. The alternative procedure, that is, modifying the geo- 
metry directly by addition of the displacement thickness, while 
quite widely used in two dimensions, becomes untenable in three 
dimensions. Primarily, this is a result of having to calculate 
and invert a new aerodynamic influence coefficient metrix at each 
iteration due to the change in geometry- Without considering the 
additional cost of smoothing the geometry and redefining the panels 
at each iteration, the additional cost of calculating and invert- 
ing the influence coefficient matrix at each iteration can be seen 
from the following illustration. 

Wing ; Kolbe and Boltz (Ref. 22 ) 

— sweep angle 45° 

— taper ratio . 5 

— aspect ratio 3 

— number of panels 510 

— computer CDC 7600 

Calculation and inversion of the aerodynamic influence coef- 
ficient matrix required 193 CP seconds on the CDC 7600. The 
iterative solution required 2.1 CP seconds. Subsequent iterations 
with the source method requried only the 2.1 CP seconds per itera- 
tion to obtain a new potential flow solution, with convergence 
acheieved in 4 iterations. The direct displacement thickness 
method would have required as a minimum an additional 193 CP 
seconds per iteration to obtain the desired solution. Assuming 
that it would also have converged in 4 iterations , an additional 
(193 X 3) = 579 CP seconds would have been needed to complete the 
viscous/inviscid calculation. On slower computers, the additonal 
time requirements would be prohibitive. 
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CALCULATION PROCEDURE 


The computer program is made up of a series of overlays, as 
shown in Figure 7. The executive program, VIP3D, controls the 
overall analysis by calling in turn the overlays containing the 
potential flow and boundary layer calculation methods. The cal- 
culation sequence is outlined as follows. 

i) The input geometry, as represented by a series of 
streamwise planar panels is lofted in Program WBPAN. 

ii) The potential flow pressvire field is computed for the 
multi-element configuration in Program WBAERO. Off- 
body pressures over flap upper surfaces are computed 
as part of the analysis. Included in the calculation 
is the force-free wake analysis. The multi-element 
wing can consist of up to four elements, a leading- 
edge slat, the main wing and double slotted flaps. 

iii) The boundary layer development is determined for each 

streamwise strip beginning inboard and proceeding span- 
wise to the wing tip. Each strip is treated as a sepa- 
rate infinite yawed wing. The laminar and turbulent 
boundary layer developments are determined for each 
element of the strip as a function of the potential 
flow pressvire distribution. Transition or laminar 
separation and turbulent separation are predicted, if 
present. Program IBL is used for the integral boundary 
layer analysis, and Program INSPAN, the finite differ- 
ence method, is used when flap upper surfaces are con- 
sidered. 

iv) Source distributions representing the boundary layer 
displacement effect are determined for each streamwise 
strip. 
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OVERLAY 



Figure 7, 


Viscous/Potential Flow Program Overlay Structure 








v) A new potential flow solution is computed taking into 
account the source distribution computed in step (iv) . 

Steps (iii) through (v) are repeated until convergence (based 
on configuration lift coefficient) is achieved. Forces and moments 
are then calculated, both for the complete configuration, and for 
the individual sections. The calculation procedure is illustrated 
in Figure 8 . 
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Figure 8. Computation Procedures 
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CALCULATIONS AND DISCUSSION OF RESULTS 


Calculations have been made to compare with experimental 
results available from a recent NASA Langley Research Center 
study. The experimental configuration consisted of a rectangular 
wing of aspect ratio = 6 mounted on a splitter plate in the wind 
tunnel. A range of angles of attack were investigated for a se- 
ries of sweep angles. Calculations made for the zero sweep case 
are shown in Figure 9 . Comparisons with experiment show some 
improvement in the predicted spanwise load distribution with vis- 
cous effects added. Predicted streamwise pressure distributions 
(fourth iteration) at three spanwise locations are in good agree- 
ment with experiment. Tip-edge vortex formation and wake roll-up 
are not modelled in this calculation, although at the angle of 
attack considered (6®) their effect would be small. There is, 
however, already some modification to the experimental pressure 
distribution at the 98.8% semispan location due to the tip-edge 
vortex at a = 6®. The boundary layer calculations were based on 
forced transition at 10% chord on the wing to correspond to a 
transition strip at 5% chord on both upper and lower surfaces of 
the wing in the experimental case. 

A second configuration investigated was the swept wing case 
of Kolbe and Boltz (Ref. . This case is of interest because 
of the low aspect ratio (3), the sweep angle (45®) and the taper 
ratio (.5). Several investigatators have used the experimental 
data to compare with their calculation methods. Of particular 
interest are the studies by Hess (Ref. ^) . His inviscid calcu- 
lations indicated a spanwise load distribution approximately 15% 
higher than experiment. He postulated that the difference was 
due to viscous effects, and attempted to prove it by, in his own 
words, adding a crude estimate of displacement thickness to the 
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Figure 9, Comparison of Measured and I 
Pressure Distributions 
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wing. The new calculations were in much closer agreement with 
experiment, and according to Hess, justified his assumption. Cal- 
culations by the present method, however, dispute this assumption. 
As shown in Figure 10, the inviscid calculation also is high by 
about 15% in spanwise loading, with only modest improvement when 
the number of spanwise strips used to model the wing is doubled. 

At the Reynolds number of the test, 18 million, the calculations 
indicate that the viscous interaction with the pressure field is 
quite small, smaller in fact than the difference between calcula- 
tions using 5 and 10 spanwise strips. Because of the high Reyivalds 
number and the symmetric, relatively thin airfoil section (7.6% 
t/c on sections parallel to the plane of symmetry) and the low 
angle of attack (a = 8®), viscous effects should be small. There- 
fore some other hypothesis is needed to explain the differences 
between theory and experiment. The explanation appears to come 
from the shape of the wing tip (rounded in the case of Kolbe and 
Bolts) and from the low aspect ratio. Experimental evidence col- 
lated by Hoerner (Ref.^ shows that rounding tip edges (lateral 
and streamwise) reduces the effective aspect ratio and the lift 
curve slope of the wing. As geometric aspect ratio decreases, the 
influence of rounding becomes more and more predominant. Since 
the inviscid calculation methods (Hess and the present method) do 
not specifically take into account the tip-edge shape, they will 
likely over-predict the lift distribution. 

A simple calculation substantiates this assumption. Sharp 
lateral-edge wings with an aspect ratio of 3 have a lift curve 
slope, dCj^/da® = .058, whereas with rounded edges, dCj^/da® = .05 

(Ref. 24). The difference in for wings having symmetric sec- 
tions at 8® angle-of-attack is approximately .065, or very close 
to the difference between theory and experiment shown in Figure 10. 
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It would appear that in order to obtain improved results for low 
aspect ratio wings, the potential flow methods must use more real- 
istic paneling in the tip region, and, perhaps even incorporate 
tip-edge vortex roll-up model. 

The multi-element capability of the program was checked using 
a wing configuration that has been tested in the Ames 40- by 80- 
foot Wind Tunnel (Ref. 2 ). The wing, consisting of the RAE 2815 
airfoil and single slotted flap cross-section had a span of 52.6 
feet, and a reference chord of 5.58 feet, giving an aspect ratio 
of 9.43. The calculations were performed using 68 chordwise 
panels and 3 streamwise strips to represent the v/ing on one side 
of the plane of symmetry. This is a rather coarse paneling, but 
it should demonstrate the trends. The calculations for lift, 
drag and pitching moment coefficients for a series of angles-of- 
attack are shown in Figure 11. The spanwise distributions of 
lift, drag and pitching moment are shown in Figure 12. The cal- 
culations which include viscous effects are determined from the 
third iteration through the program. A second calculation was 
made with the wing and flap wakes allowed to relax in a plane 
parallel to the plane of symmetry. This resulted in a slightly 
reduced lift coefficient as indicated in Figure 11. Further work 
will include detailed comparisons between the predictions of this 
theoretical method and the data from the Ames experimental program. 
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WING LIFT COEFFICIENT 



Figure 11, Predicted Lift, Drag and Pitching 
Wing with Single Slotted Flap 
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Figure 12. Predicted Spanwise Distributions of Lift, Drag and Pitching Moment 
for an R/.E 2 815 Wing with Single Slotted Flap. 



CONCLUSIONS AND RECOMI4ENDATIONS 


The three-dimensional viscous/potential flow interaction 
analysis program described in this report is currently the only 
one of its kind. Its present capability includes the analysis of 
multi-element wings having full span slats and slotted flaps. A 
:force-free (relaxed) wake model is included in order to improve 
the prediction of surface pressures on the individual wing com- 
ponents. Wake roll-up and tip-edge vortex roll-up are not yet 
modeled, but will be required in the future if partial span slats 
or flaps are to be considered. 

The method in its present form has been compared with only a 
very small data set; however, the various programs which make up 
the complete program have been used independently or in combina- 
tion with other programs for some time. In general, the methods 
are stable and have demonstrated good agreement with experiment 
over a wide range of applications. It is concluded that the pre- 
sent program will also prove to be very useful over a wide range 
of three-dimensional wing applications. 

A specific conclusion can be drawn in regards to the use of 
sources to represent the displacement effect of the boundary layer 
on the potential flow. In three dimensions, the computational 
superiority of this approach is clearly what makes such a calcula- 
tion procedure practical. The direct addition of displacement 
thickness with the accompanying necessity to distribute and smooth 
the new surface, and svibsequently reinvert the aerodynamic influ- 
ence matrix at each iteration is almost prohibitive of computer 
time . 

It is recommended that once the close- approach problem of a 
vortex wake interacting with the vortex elements of a downstream 
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surface is solved, the calculation procedure be expanded to inclvide 
models for complete wake roll-up and tip-edge vortex roll-up. 
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APPENDIX I 


PANEL GEOMETRY CALCULATION PROCEDURE 

The analytical procedure presented here follows closely the 
method first developed in Reference 14 . A quadrilateral surface 
element is described by four corner points, not nesessarily lying 
in the same plane, as shown in the sketch. Note that the number- 
ing convention of the corner points differs from that used in the 
preceding text. The quadrilateral element is approximated by a 
planar panel as follows: 
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The coordinates in the reference coordinate system are identified 


by their subscripts. 

The components 

of the diagonal vectors T^ 

and T are 
2 
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We may now obtain a vector N (and its components) by taking the 
cross product of the diagonal vectors. 


-j- -> 

N = Tj X Tj 


N = T, T - T T 

X 2y 12 ly 22 


N = T T - T T 
y IX 22 2 X 12 


N = T T 
2 2X ly 


T T 
IX 2y 


The unit normal vector, n, to the plane of the element is taken as 

-V 

N divided by its own length, N (direction cosines of outward unit 
normal) . 
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where N 


= [ 


1 ^ 

X y z J 


The plane of the element is now completely determined if a point 
in this plane is specified. This point is taken as the point 
whose coordinates, x, y, z, are the averages of the coordinates of 
the four input points. 
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Now the input points will be projected into the plane of the ele- 
ment along the normal vector. The resulting points are the corner 
points of the quadrilateral element. The input points are equi- 
distant from the plane, and this distance is 


d . = 

! 

(X - 

X ) + ny{y 

“ y ) +‘ n^ (z - z ) 1 

The coordinates 

of 

the 

corner points 

in the reference coordinate 

system are given by 




II 

^k 

+ 

(-l)^'^^n d 

X 



^k 

+ 

(-l)^'^^nyd 

k = 1, 2,' 3, 4 


^k 

+ 

(-l)^'^^n^d 



The element coordinate system is now constructed. This re- 
quires the components of three mutually perpendicular unit vectors, 
one of which points along each of the coordinate axes of the sys- 
tem, and also the coordinates of the origin of the coordinate sys- 
tem. All these quantities must be given in terms of the reference 
coordinate system. The unit normal vector is taken as one of the 
xinit vectors, so two perpendicular unit vectors in the plane of 
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the element are needed. Denote these unit vectors t and t . The 
■+ -»- 12 

vector tj is taken as divided by its own length, T^, i.e.. 
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The vector t 

1 

is 

defined 

by t = n X t , so that 

its components 

are 
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The vector t^ 

is 

the unit vector parallel to the 

X or C axis of 


the element coordinate system, while t^ is parallel to the y or n 

axis, and n Is parallel to the z or C axis of this coordinate sys- 
tem. 


The corner points are now transformed into the element coordi- 
nate system based on the average point as origin. These points 
have coordinates y^^, in the reference coordinate system. 
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Their coordinates in the element coordinate system with this 
origin are denoted by 0. Becuase they lie in the plane of 

the element, they have a zero z or 5 coordinate in the element co- 
ordinate system. In the (5,ri) coordinate system, the corner points 
of the element are: 

5k ' - =''k> ‘.y<? - yp - ^k> 

\ - ='k> + '^2y<y - ^k' 

These corner points are taken as the corners of a plane quadrila- 
teral . The origin of the element coordinate system is now trans- 
ferred to the centroid of the area of the quadrilateral. With the 
average point as origin, the coordinates of the centroid in the 
element system are: 

^ n, - n, -V 5,(n. - n.)] 

s - - i 

These are subtracted from the coordinates of the comer points in 
the element coordinate system based on the average point as origin 
to obtain the coordinates of the corner points in the element co- 
ordinate system based on the centroid as origin. Accordingly, 
these latter coordinates are 

k = 1, 2, 3, 4 
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Since the centroid is to be used as the control point of the ele- 
ment, its coordinates in the reference coordinate system are re- 
quried. These coordinates are: 




y = y + tE + tn 
■^0 ^ ly^o ay'o 


z = 2 + t E + t,n 

0 1 0 2 z ' 0 


Finally, the area of the quadrilateral is: 

A = ^(C -5){n - n) 

Z 3 12 4 


One additional rotation of the C, n axis is performed in 
order to ensure that the 5, ? phase is parallel to the reference 
X axis of the configuration. 

The direction cosines of the axial vector, A, and transverse 
vector, m, are given in terms of the normal vector, n, as follows 


n 

I = — 

X m 


A = 0 

y 


m 


n i 

y z 


m 


= V 



m. 


-n S, 

y X 
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APPENDIX I 


INFLUENCE OF TIP AND TRAILING-EDGE VORTICES 


The perturbation velocity components induced by the edge and 
trailing vortex systems associated with a vortex panel are derived 
in Reference Formulas for constant and linearly varying vor- 

tices are listed below. The perturbation velocities are given at 
an arbitrary point, P(x,y), and induced by a vortex which lies 
along the x axis between x = 0 and x = “. 

Constant Line Vortex 


u = 0 



' ^ * S J 


where 


= y 


and 

d^ = x^ + 
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Linearly Varying Line Vortex 


u = 0 



For vortex panels having a constant distribution of vorticity, 
and non-parallel leading and trailing edges, a constant vortex 
sheet is generated and contained between the two trailing vortices 
in the wake. The perturbation velocities induced by this vortex 
sheet, which originates along the panel trailing edge and extends 
downstream to infinity are given belov/: 


where 


u = 
v = 

0 



w = 



«.a> 

3 “ 

x - 

log 

c, + V (X - 

2 

X - 

C3 + V (X - C3)" 

z 

+ ^3 



and j / j ' defined following equations ( 3 ) in the 

in the main text. The subscripts 3 and *♦ refer to the outboard 
and inboard trailing edge corners of the vortex panel. 
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All of the above formulas are v/ritten assuming that the wake 
lies in the plane of the vortex panel, and that the edge and trail- 
ing vortices are parallel to the ^ axis of the panel. If the panel 
side edges or trailing vortices lie in the panel ? ,ri plane, but 
are inclined at an angle, 6 » with respect to the ? axis, the follow- 
ing transformation applied: 


u 

= u' 

cos 6 - 

v' sin 6 

v 

= v' 

cos 6 + 

u* sin <5 

w 

= w* 



v' 

and w 

' are the 

three components of velocity calcula 


ted in a primed coordinate system aligned with the vortex under 
consideration. 
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